Genomic divergence, local adaptation, and complex demographic history may inform management of a popular sportfish species complex

Abstract The Neosho Bass (Micropterus velox), a former subspecies of the keystone top‐predator and globally popular Smallmouth Bass (M. dolomieu), is endemic and narrowly restricted to small, clear streams of the Arkansas River Basin in the Central Interior Highlands (CIH) ecoregion, USA. Previous studies have detected some morphological, genetic, and genomic differentiation between the Neosho and Smallmouth Basses; however, the extent of neutral and adaptive divergence and patterns of intraspecific diversity are poorly understood. Furthermore, lineage diversification has likely been impacted by gene flow in some Neosho populations, which may be due to a combination of natural biogeographic processes and anthropogenic introductions. We assessed: (1) lineage divergence, (2) local directional selection (adaptive divergence), and (3) demographic history among Smallmouth Bass populations in the CIH using population genomic analyses of 50,828 single‐nucleotide polymorphisms (SNPs) obtained through ddRAD‐seq. Neosho and Smallmouth Bass formed monophyletic clades with 100% bootstrap support. We identified two major lineages within each species. We discovered six Neosho Bass populations (two nonadmixed and four admixed) and three nonadmixed Smallmouth Bass populations. We detected 29 SNPs putatively under directional selection in the Neosho range, suggesting populations may be locally adapted. Two populations were admixed via recent asymmetric secondary contact, perhaps after anthropogenic introduction. Two other populations were likely admixed via combinations of ancient and recent processes. These species comprise independently evolving lineages, some having experienced historical and natural admixture. These results may be critical for management of Neosho Bass as a distinct species and may aid in the conservation of other species with complex biogeographic histories.


| INTRODUC TI ON
The convention of classifying organisms into discrete taxonomic units, typically "species," before they can receive conservation priority (Beheregaray & Caccone, 2007;Isaac et al., 2004) is often politically charged and ignores the biological reality that any arbitrary unit is composed of nested genetic groups. Reciprocally monophyletic lineages at the highest tier of differentiation are made up of metapopulations; metapopulations may consist of many populations; populations are divided into subpopulations and finally into pedigrees. Each of these levels may contain valuable allelic polymorphisms (Lawson, 2013;Praebel et al., 2013) which, in quickly fluctuating environments, could provide the raw material for cladogenesis (Hendry, 2017). Delineating intraspecific variation at the genomic level may therefore aid in biodiversity conservation, especially when that variation is cryptic and perhaps overlooked due to behavioral traits or convergent morphology (Culver et al., 1995;Culver et al., 2009;Schluter, 1996). It is equally crucial to ascertain the eco-evolutionary context leading to contemporary diversity to predict how species will adapt in the changing world.
Characterizing amounts and patterns of genetic diversity (e.g., allelic richness, allele frequency differentiation, and admixture) and their underlying causal mechanisms (e.g., selection, drift, and gene flow) is challenging for freshwater riverine wildlife. The onedimensional, dendritic configuration and abiotic heterogeneity of stream ecosystems, including variable flow rates, depths, temperature gradients, nutrient concentrations, and allochthonous and autochthonous inputs (Barthel et al., 2008;Lytle & Poff, 2004;Vannote et al., 1980), create diverse conditions and restrict movement to relatively narrow corridors within watersheds, setting the stage for population structure (Herdegen et al., 2014;Jacobsen & Hansen, 2005;Puebla, 2009;Ward et al., 1994), local adaptation, and potentially distinct demographic histories among populations.
Life history and behavior, such as habitat use, dispersal, and reproduction, may also contribute to eco-evolutionary dynamics.
Fish species that are valued for angling or aquaculture face exceptionally complex environmental pressures because they may also be subjected to human-mediated introductions (Hohenlohe et al., 2013).

Phylogeography of North American endemic black basses
(Micropterus) is only partially understood but is likely to have been shaped by both natural and anthropogenic forces. One of the most economically important and globally popular black bass species, the Smallmouth Bass (M. dolomieu), occupies a native range extending from the Laurentian Great Lakes in southeastern Canada to the Central Interior Highlands, USA (CIH). Such a broad, ecologically variable distribution and high levels of range-wide diversity (e.g., Borden & Krebs, 2009) have made it especially difficult to resolve the species' taxonomy. Biologists historically recognized two subspecies, with one encompassing the central and eastern portion of the range (Northern Smallmouth Bass, M. d. dolomieu;Hubbs & Bailey, 1940) and another being restricted to the Arkansas River Basin (Neosho Smallmouth Bass, M. d. velox;Hubbs & Bailey, 1940).
Lack of genome-wide assessments prevented fine-scale resolution of differentiation between the subspecies and precluded their designation as independently evolving lineages.
A recent phylogenomic study of the black bass genus provided compelling evidence of genomic divergence between the subspecies (Kim et al., 2022). The authors ultimately elevated the Neosho Smallmouth Bass to species rank (Neosho Bass; M. velox), consolidating the Northern Smallmouth Bass as synonymous with Smallmouth Bass (M. dolomieu). Recent investigations of morphological and ecological divergence have largely affirmed this taxonomic revision Gunn et al., 2020;Hubbs & Bailey, 1940;Miller & Brewer, 2021. However, other genetic studies have revealed considerable population structure within both native ranges (Gunn et al., 2020;Long et al., 2021;Stark & Echelle, 1998;Taylor et al., 2018) and substantial, heterogeneous admixture within the Neosho range (Gunn et al., 2020;Taylor et al., 2018), in contrast to the previously assumed scenario of two diverging allopatric lineages.
These studies suggest a more complex dichotomy of differentiation and gene flow which varies across clades, streams, and populations, making it challenging to discern evolutionarily significant units (Moritz, 1994) and proceed with conservation and protection priorities.
The most recent time-calibrated phylogeny of the black basses Kim et al., 2022) dates the split of Smallmouth Bass from its sister clade to between 4 and 6 million years ago.
Advancement of the glacial front to the last glacial maximum (~22-18 thousand years ago), which coincides with the parapatric convergence of the Smallmouth Bass and Neosho Bass ranges, would have pushed fish into southern refugia (Borden & Krebs, 2009). Later recession may have altered the topography enough to sever river connections, creating opportunities for vicariant speciation. Recession may have also joined rivers through erosion, allowing for dispersal and subsequent gene flow (Berendzen et al., 2008;Near et al., 2001;Near & Keck, 2005;Ray et al., 2006;Zink, 1997). This aligns with the fact that the CIH is an endemism hotspot (Soltis et al., 2006) for freshwater fishes (Cross et al., 1986;Lundberg et al., 2010;McAllister et al., 1986;Robison, 1986) and lends anecdotal support

T A X O N O M Y C L A S S I F I C A T I O N
Applied ecology, Biogeography, Conservation genetics, Ecological genetics, Evolutionary ecology, Genomics, Population genetics, Taxonomy to the possibility of greater inter-and intraspecific diversity in black basses in this ecoregion.
In addition to historical geological and ecological shifts, post-Pleistocene genetic structure may be influenced by contemporary processes. Smallmouth and Neosho Bass populations exhibit inconsistent dispersal and migratory behavior. Some are sedentary (Funk, 1957) due to philopatry and nest-site fidelity (Ridgway et al., 1991), while others are seasonally potamodromous (Funk, 1957;Gowan et al., 1994;Lyons & Kanehl, 2002); these behaviors may to some degree vary by species (Miller & Brewer, 2021).
Smallmouth and Neosho Bass are also extremely economically valuable; hatchery-raised Smallmouth Bass have been introduced around the globe for recreation, trophy angling, and as a food source (Brewer & Orth, 2015;Iguchi et al., 2004;Robbins & MacCrimmon, 1974;Stark & Echelle, 1998). A genetic strain of Smallmouth Bass derived from the Cumberland River drainage was used to stock the Illinois River system within the Neosho Bass native range in the early 1990s . However, this single known stocking event does not adequately explain signatures of substantial, heterogeneous admixture in Neosho streams (Gunn et al., 2020). Unreported or inadvertent introductions may be responsible; otherwise, admixture may be a natural byproduct of stream piracy (Branson, 1963) or recent flooding. Both scenarios substantially impact the genetic integrity of the species in this region.
Elevating the Neosho Bass to species status has profound implications for conservation and management of economically and ecologically valuable populations in a popular sportfish species complex in the CIH. It is therefore critical to understand the extent of divergence, the diversifying mechanisms generating inter-and intraspecific diversity, and the homogenizing forces potentially eroding adaptive variation to inform effective strategies for long-term viability. Genomic sequencing technologies provide the resolution and power to study highly structured populations in complex physical environments which may be susceptible to gene flow. We harnessed reduced representation sequencing (ddRAD-seq; Peterson et al., 2012) to resolve genomic diversity, local adaptation, and demographic history in the Smallmouth Bass and Neosho Bass, representing some of the world's most popular game fisheries. We specifically examined: (1) phylogenetic hypotheses between and within species; (2) differentiation between and within species at outlier loci; and (3) alternative admixture scenarios using a model-testing framework, in which we inferred joint demographic histories from population-specific site frequency spectra. We expected that the Smallmouth and Neosho Bass would be reciprocally monophyletic and that populations would be nested within species. We expected that outlier loci would more strongly differentiate species than populations. Finally, we expected that the genetic architecture of all admixed populations within species would be best explained by very recent gene flow, implicating anthropogenic introductions.

| Sample collection and genomic DNA preparation
We obtained 95 black bass samples, representing the Neosho Bass (N = 66) from 13 streams throughout the Arkansas River Basin (ARB), and the Smallmouth Bass (N = 25) from three tributaries of the White River (WRT), two tributaries of the Missouri River (MRT), and Skiatook Lake (LAKE), an impoundment in northeastern Oklahoma situated outside the native range of Smallmouth Bass that was stocked with a hatchery-reared strain colloquially known as "Tennessee lake-strain" sourced from the Cumberland River drainage (CIH; Figure 1a and Table 1; Table S1; Gunn et al., 2020). For phylogenomic comparison, we used four Spotted Bass (M. punctulatus) from the ARB as an outgroup ( Table 1). Including Spotted Bass, we sampled from 20 total stream or impoundment sites.
We extracted high molecular weight DNA (gDNA) from ~25 mg fin clips excised from the upper caudal fin using the DNeasy Blood and Tissue kit (QIAGEN, Germantown, MD). Fin clips were coarsely chopped with sterile razor blades, digested with proteinase K, and treated with 4 μl RNase-A. Extracts were diluted to ~20 ng/μl with ddH 2 O, arranged on a 96-well plate in 50 μl (~100 μg gDNA) aliquots and stored at −20°C before library preparation. We included a single, no-DNA negative control sample containing only ddH 2 O.

| Library preparation and sequencing
Library preparation and sequencing for ddRAD-seq were completed at Floragenex, Inc. (Eugene, OR) according to a modified Sequencebased Genotyping protocol outlined in (Truong et al., 2012).
Approximately 500 ng genomic DNA was digested with PstI and MseI at 65°C for 1 h, followed by ligation of paired-end P5 PstI and AFLP MseI adaptors at 37°C for 3 h. Unique, 5-base pair (bp) barcodes were included in the P5 PstI adaptor for individual sample identification. PCR was performed using the parameters listed in (Truong et al., 2012). Samples (N = 95) were pooled and sequenced with 1 × 95 bp chemistry on a single lane of the Illumina HiSeq 4000.

| Bioinformatic processing
Sequence filtering, clustering, alignment, and assembly were completed by Floragenex, Inc. (Eugene, OR) according to the RADseq processing and variant detection pipeline (Lozier, 2014). Clusters from the sample with the highest number of unique clusters (AR21; Table S1) were used to assemble a de novo reference sequence, which was aligned back to itself to minimize paralogs and to which all Smallmouth and Spotted Bass individuals were aligned (Lozier, 2014).

| SNP discovery and filtering
Clusters were processed into RAD tags (95 bp sequences), and single-nucleotide polymorphisms (SNPs) were called using SamtoolS v.0.1.16 (Li et al., 2009)  Variants were initially filtered based on individual read depth; only sequences with a minimum of 15X coverage were retained. We removed samples from the dataset that had greater than 20% missing genotype calls across SNPs. We then removed SNPs with phred quality scores less than 20 (Liao et al., 2017) and greater than 20% missing genotype calls across samples (e.g., Lavretsky et al., 2019).
We then generated two datasets, one in which RAD tags were thinned to retain a single SNP (to reduce the likelihood of linkage between variants in phylogenomic analyses), and one in which RAD tags were not thinned (for fine-scale population delimitation requiring multiple SNPs per RAD tag to increase computation power).
For both the thinned and nonthinned datasets, we removed remaining SNPs with minor allele count of two or less (equivalent to a minor allele frequency of ~0.011). We created SNP tables using Gatk we eliminated any remaining paralogs by omitting SNPs that were heterozygous in greater than 45% of samples.

| Lineage divergence
We screened the full dataset for individuals of putative admixed ancestry, i.e., those of Smallmouth Bass × Spotted Bass or Neosho Bass × Smallmouth Bass hybrid origin, to limit gene flow influence on the assessment of lineage diversification. Our full VCF file was F I G U R E 1 Species geographic ranges, sampling sites, and distinct evolutionary lineages. (a) Native ranges of the Smallmouth Bass (Micropterus dolomieu; light gray) and the Neosho Bass (M. velox; dark gray), with representative illustrations. (b) Maximum-likelihood phylogeny for putatively pure (p-Pure) Spotted Bass, Smallmouth Bass, and Neosho Bass, with black and gray boxes at nodes indicating 100% and > 90% bootstrap support, respectively, and population structure results for K = 3, K = 4, and K = 5, with major lineages and sample sites labeled corresponding to individual samples. (c) 10-fold crossvalidation error results for admixture analysis. (d) Sampling sites (numbered as in Table 1) within the Central Interior Highlands (CIH) for Smallmouth and Neosho Bass colored by distinct evolutionary lineages. Sites of putative admixed origin based on preliminary admixture analysis (p-Admixed) are indicated as white circles; empty white circles indicate sites where all individuals were of putatively admixed origin, and stars indicate sites where nearly all individuals were of putatively admixed origin.
converted to binary pedigree (BED) format in a high-contig build of Plink v.1.90 (Chang et al., 2015). We estimated ancestry proportions (q) for individual fish in the program admixture v.1.3.0 (Alexander et al., 2009), inferring the optimal number of K clusters by minimizing 10-fold cross-validation error for K = 1-20 (number of stream sites plus one additional cluster for the Spotted Bass outgroup). We used stringent criteria to determine pure or admixed origin: individuals were considered putatively pure if q ≥ 0.95 for one inferred cluster at the optimal K. (Thongda et al., 2020). While minor cluster membership of 0 ≤ q ≤ 0.05 for an individual may indicate historical introgression, we retained individuals in this range for phylogenomic comparison to avoid excessively limiting sample sizes. Other studies have resolved lineage divergence using individuals with minor ancestry of 0 ≤ q ≤ 0.25 . Hybrids of Smallmouth Bass or Neosho Bass with Spotted Bass were removed from downstream analyses. The dataset was separated into two subsets: (1) putatively pure individuals ("p-Pure") and (2) putatively admixed individuals ("p-Admixed").
We used the p-Pure dataset to investigate phylogenomic relationships among and within Smallmouth and Neosho Bass. We first assessed allele frequency differentiation in the population structure program admixture, choosing K by minimizing 10-fold cross-validation error for K = 1-20. We conducted a parallel phylogenomic analysis using maximum-likelihood methods in the SnPhylo pipeline (Edgar, 2004;Felsenstein, 1989;Lee et al., 2014;Schliep, 2011;Zheng et al., 2012). We ran 10,000 bootstrap replicates (−b) using Spotted Bass as an outgroup (−o) and setting a linkage disequilibrium threshold (−l) of 0.1. A consensus tree was constructed in fiGtree v.1.4.2 (Rambaut, 2010) and aligned with results from admixture.

| Population discovery
To delimit populations and assess connectivity among stream sites, we examined haplotype similarity among Smallmouth and Neosho Bass individuals in fineradStructure v.0.3.2 (Malinsky et al., 2018).
We used our full, nonthinned SNP dataset and concatenated SNPs on the same RAD tag to form haplotypes in the radPainter package (Malinsky et al., 2018). We calculated a co-ancestry matrix with 100,000 burn-in steps (−x), 100,000 Markov chain Monte Carlo (MCMC) iterations (−y), and thinning (−z) every 1000 iterations for the full sample set, including p-Pure and p-Admixed individuals. We also calculated co-ancestry matrices separately for the p-Pure and p-Admixed groups to eliminate bias due to multiple ancestry; results one Tennessee Lake-strain-stocked lake population (LAKE) before data filtering.
from the separate matrices were used for population delimitation.
We used a full hill-climbing tree-building method to construct trees, running 10,000 iterations (−x), providing no value for the initialization parameter (−T). Individuals were collapsed into populations if they formed blocks of high co-ancestry along the diagonal of the co-ancestry matrix and if they were monophyletic at deeper nodes in the tree.

| Adaptive divergence
We explored the role of adaptive divergence influencing variation among and within species by scanning for outlier SNPs putatively showing high (diversifying selection) or low (balancing selection) differentiation among populations (p-Pure and p-Admixed). Outlier analyses are effective in identifying SNPs that deviate from null expectations of allele frequencies under an island model of migration (Lewontin & Krakauer, 1973). However, they are often prone to high rates of false positives, especially when the studied populations are distributed spatially in one-dimensional, steppingstone arrangements, as is the case for riverine fish species (Bierne et al., 2013;Fourcade et al., 2013). We alleviated bias from potential false-positive results (Excoffier et al., 2009;Jakobsson et al., 2013;Jost, 2008;Nei & Maruyama, 1975;Robertson, 1975) by comparing outliers from different statistical analyses.
We conducted genome scans hierarchically at the black bass species level, among species in the Smallmouth Bass species complex, and among populations within species (Chen et al., 2016) to reduce the effect of population structure on outlier detection. We used our full, thinned SNP dataset, dividing individuals into four groups: (1)  For each hierarchical analysis, we used default MCMC parameters in BayeScan, retaining only SNPs with logged posterior probability greater than 1.5, deemed "very strong" support for selection (Foll & Gaggiotti, 2008). For PCAdapt, we tested K = 1-20 principal components (PCs). To determine the optimal number of PCs, we assessed Scree plots and selected the number of PCs based on Cattell's Rule (Luu et al., 2017). We generated p-values for all SNPs, applying a Bonferroni correction for multiple tests.
Final sets of outlier SNPs were created by merging candidate outliers from BayeScan and PCAdapt. To assess neutral differentiation (drift), we also created datasets with only shared neutral SNPs (nonoutliers).
We plotted samples according to our a priori population designations (Miller et al., 2020)

| Admixture mapping
Admixture signatures could be due to gene flow or incomplete lineage sorting. We tested for evidence of these processes using mixmaPPer v.2.0 (Lipson et al., 2013(Lipson et al., , 2014. Populations with shared alleles due to incomplete lineage sorting are inferred as nonadmixed, whereas those with a history of admixture postdivergence are inferred as admixed. We used the full dataset with all inferred populations and created input files in eiGenSoft v.7.2.1 Price et al., 2006). mixmaPPer uses physical and genetic linkage to calculate genetic drift after admixture. Since our data were not mapped to a reference genome, we did not have linkage information and therefore did not infer precise divergence and mixing times using the drift units generated. To identify significantly admixed populations and fit them to a scaffold tree, we assumed independence (no genetic or physical linkage) of all SNPs (given that we filtered for one SNP per RAD tag during bioinformatic processing) and therefore generated arbitrarily large physical and genetic distance values for each SNP according to the custom formula: where d = physical/genetic distance, x = numerical label of the RAD tag (1, 2, …), and y = nucleotide coordinate of SNPs within RAD tags.
Moment statistics were calculated using 1000 bootstrap replicates over 50 cM blocks, and the scaffold tree was constructed using 10,000 data subsets. For populations not included in the scaffold tree, we tested the fit of two-way admixtures between all pairs of nonadmixed parents (sources), running 100 bootstrap replicates. Significantly admixed populations were used in demographic analyses.

| Demographic history
We explored potential demographic scenarios driving observed admixture patterns between p-Pure and p-Admixed populations, testing nine two-population diversification-based demographic models (Portik et al., 2017)  two-population models allow for divergence with and/or without migration between focal populations. Candidate scenarios differed in the timing of migration, i.e., ancient or due to recent secondary contact, and directionality of migration, i.e., symmetric or asymmetric. Model descriptions are given in Table S2, and parameters estimated are described in Table S3.
Demographic inference in a i assumes SNPs are unlinked and neutral (Gutenkunst et al., 2009). Thus, we used only neutral SNPs ascertained from adaptive divergence analysis, converting SNP data for each population pair into folded 2D joint site frequency spectra It is theoretically feasible to convert a i parameter estimates to measures of migration rates, divergence times, and population sizes.
However, since we do not know the true demographic history of these populations, it is possible none of our candidate models fully explain genetic diversity. Additionally, parameter conversion should be conducted with a bootstrapping procedure to quantify uncertainty (Gutenkunst et al., 2009;Portik et al., 2017), and we have rel- We removed three samples (GRSPB23, ER05, and BFORK32; Table S1) from the dataset that had greater than 20% missing genotype calls ( Figure S1). After all filtering, the final dataset (N = 92, 24 Smallmouth Bass, 64 Neosho Bass, and 4 Spotted Bass) contained 50,828 SNPs for downstream analyses.

| Lineage divergence
Population structure results for all samples were supported at K = 4 by 10-fold cross-validation (CV error = 0.253; Figure S3), revealing that all Smallmouth Bass were putatively of pure origin, but 64% of all Neosho Bass (N = 41) were admixed while 36% (N = 23) were putatively pure ( Figure S4). One Neosho Bass was identified as a likely Spotted Bass hybrid (BFC10 ; Table S1) and was removed from downstream analyses. One or more p-Admixed individuals were identified in all but three Neosho Bass sample sites (Sites 8, 11, and 17). In six Neosho sampling sites (Sites 7,9,10,12,18,and 19), all individuals were putatively admixed ( Figure S4).

| Population discovery
Co-ancestry analysis of the full sample set did not resolve populations corresponding to rivers, instead showing paraphyly among individuals collected from the same site ( Figure S5)

| Adaptive divergence
Scanning our SNPs with BayeScan across all black bass samples re- Neosho Bass only, we found 32 candidate diversifying SNPs and 50,796 neutral SNPs ( Figure S6c). Among Smallmouth Bass only, only six candidate SNPs were found to be under substantial selection, and we found 50,822 neutral SNPs ( Figure S6d).
In PCAdapt, we retained two PCs for analysis with all black bass samples ( Figure S7a

| Admixture mapping
All p-Pure populations (WHITE, MISS, SKIA, MIDARK, and LMULB) were found to be nonadmixed and formed the branch tips of the admixture tree ( Figure 4). All p-Admixed populations (ELK, ILLI, UPPARK, and BAYOU) were significantly admixed based on f 3 statistics (Table S4). All admixed populations were inferred to be parented by the MIDARK population of Neosho Bass. The ELK, BAYOU, and UPPARK populations were all admixed with the WHITE population of Smallmouth Bass lineage ( Table 2 and Figure 4). The ILLI population was admixed with the SKIA population (Table 2 and Figure 4). In each case, sources were inferred with 100% bootstrap support ( Table 2).

| Demographic history
All nine demographic models tested in a i were successfully optimized for the ELK and WHITE (Figure S10a   Table 3 and Figure 5a).
The best-fitting model for ILLI and SKIA was SCAM (Table S2).
Migration rate from SKIA to ILLI (m 12 = 3.434) was higher than in the opposite direction (ILLI to SKIA; m 21 = 0.276). The estimate for the period of isolation following divergence (T 1 = 1.143) was substantially greater than the time since secondary contact (T 2 = 0.091; Table 3; Figure 5b).
The AM2E model (Table S2) was the most suitable for UPPARK and WHITE, ( Table 3; Table S5; Figure 5c). In both the first and second epochs, migration rate from WHITE to UPPARK (m 12a = 1.656 and m 12b = 1.413, respectively) was greater than the opposite direction (m 21a = 0.096 and m 21b = 0.011). Timing of each epoch with respect to species divergence was variable, although the second epoch (T 2 = 16.128) was estimated to last longer than the first (T 1 = 3.076; Table 3 and Figure 5c).
The most suitable model for BAYOU and WHITE was SCAM (Table S2). Migration rate from WHITE to BAYOU (m 12 = 19.782) was much higher than in the opposite direction (BAYOU to WHITE; m 21 = 0.775). The estimate for the period of isolation following divergence (T 1 = 0.750) was substantially greater than the time since secondary contact (T 2 = 0.019; Table 3 and Figure 5d).

F I G U R E 2
Co-ancestry and phylogenomic relationships between Smallmouth Bass (Micropterus dolomieu) and Neosho Bass (M. velox) inferred in fineradStructure for (a) p-Pure samples, and (b) p-Admixed samples. We used the full non-thinned dataset (98,659 SNPs) to generate SNP haplotypes for co-ancestry and phylogenomic assessment. Colors within the co-ancestry matrices reflect the extent of coancestry between adjacent samples, with white and blue representing low co-ancestry, and pink and black representing high co-ancestry. Residual estimates and distributions for all best-fitting models are given in Figure S14.

| DISCUSS ION
Recent studies have attempted a more holistic, multidimensional approach to population genomic investigations which seek not only to reveal the scope of diversity, but also to identify and quantify the strength of evolutionary forces acting on groups of organisms (Bangs et al., 2020;Ebersbach et al., 2020;Portik et al., 2017). By applying this approach and accounting for the many spatially and temporally dynamic processes responsible for phylogeographic patterns, we detected and described the complex genomic diversity of one of the world's most ecologically and economically valuable freshwater sportfish species.
Our phylogenomic and population genomic analyses showed that the Smallmouth and Neosho Bass are highly diverged, forming  (Gunn et al., 2020). This finding is surprising but makes sense given that endemism in the White River has been observed in other fishes (Roe et al., 2008). A second lineage within the Smallmouth Bass inhabits tributaries of the Missouri River but also includes a reciprocally monophyletic Tennessee lake-strain population in Skiatook Lake, Oklahoma, USA. Although Skiatook Lake is known to have been stocked with the Tennessee lake-strain , our data confirm it is derived from Smallmouth Bass.
Within the Neosho range, we found a nonadmixed, monophy- Concordant with their independent evolutionary trajectories, the Smallmouth and Neosho Bass are highly differentiated at multiple outlier SNPs across the genome, which could potentially be explained by directional selection in their local environments.
Previously described morphological differences between the species, including head length and number of soft dorsal fin rays (Gunn et al., 2020;Hubbs & Bailey, 1940), have provided support for adaptive diversity. However, this is the first evidence to explicitly show a genomic basis for local selection, and potentially local adaptation, rather than phenotypic plasticity. Additionally, two populations in the Neosho range, one in the Illinois Bayou River and Big Piney Creek and another in Lee Creek and the Mulberry River, were strongly diverged from all other Neosho populations and were differentiated F I G U R E 4 Admixture drift tree for Smallmouth Bass (Micropterus dolomieu) and Neosho Bass (M. velox) populations inferred in mixmaPPer. Black labels plotted on the tips of the scaffold tree with corresponding-colored squares represent pure, unadmixed populations (as determined by a two-way f 3 test). Colored and shadowed branches and labels mapped onto the scaffold tree represent significantly admixed populations originating from their respective parents. The scale for branch lengths is in drift units (D) in which D is roughly equal to 2F ST .  (Kirkpatrick & Barton, 2006), genetic hitchhiking due to selective sweeps (e.g., Kim & Maruki, 2011), purging of deleterious alleles (Pannell & Charlesworth, 2000), demographic history (Lotterhos & Whitlock, 2014), and heterogeneity in recombination rates (Roesti et al., 2012). Fine-scale genome mapping of candidate selected loci should be conducted to disentangle these effects before local adaptation can be definitively inferred.
Phenotypic differentiation may occur along ecological clines (Conover et al., 2009;Savolainen et al., 2013), and clines can be especially steep among or within rivers (Schlosser, 1991;Vannote et al., 1980). There can be substantial thermal heterogeneity or variation in hydrological factors such as flow rate, depth, and frequency and magnitude of stochastic disturbances, that is, flooding and drought (Barthel et al., 2008;Lytle & Poff, 2004). Thus, fish in different populations may be adapted to specific combinations of variables (Aitken et al., 2008;Davis & Shaw, 2001;Franks & Hoffmann, 2012), and, at the genomic level, such adaptations may be highly polygenic. Both the Illinois Bayou River and Big Piney Creek flow southward through the Boston Mountains of northern Arkansas, USA, before emptying into the Arkansas River. These streams are warm and flow along steep topographical gradients; water temperatures may exceed 30°C (Hafs et al., 2010) in summer, approaching the critical swimming maximum temperature (35°C) for fry (Larimore & Duever, 1968). Although many streams in the Neosho range flow continuously (Robison & Buchanan, 1984), Boston Mountains streams are not spring-fed and partially dry during summer and autumn months (Hines, 1975), leaving only deep, isolated pools for refuge (Hafs et al., 2010). Neosho Bass in the Boston Mountains may therefore experience occasional isolation and may be well-adapted to extreme temperature and flow regimes, warranting further investigation and conservation actions given climate projection for warmer temperatures and increased drought intensity for the region (Sharma & Jackson, 2008  . The thermal tolerance of Smallmouth Bass occupying the Boston Mountains is unknown. We identified many SNPs that were selectively neutral among Smallmouth and Neosho Bass at both the species and population levels, indicating drift. Some Smallmouth Bass individuals are sedentary (Funk, 1957) as a result of philopatry (Ridgway et al., 1991), which may contribute to reproductive isolation between populations and allow for random fixation of distinct alleles. While some fish may also be migratory, when migration occurs, it is typically seasonal (Funk, 1957;Gowan et al., 1994;Lyons & Kanehl, 2002). Abbreviations: θ, genetic diversity considering only polymorphic SNPs; AM, divergence with asymmetric migration; AM2E, divergence, asymmetric migration between two distinct epochs; nu 1 , size of population 1; nu 2 , size of population 2; m 12 , continuous migration rate from population 2 to population 1; m 12a , migration rate from population 2 to population 1 during first epoch; m 12b , migration rate from population 2 to population 1 during second epoch; m 21 , continuous migration rate from population 1 to population 2; m 21a , migration rate from population 1 to population 2 during first epoch; m 21b , migration rate from population 1 to population 2 during second epoch; SCAM, divergence, isolation, and secondary contact with asymmetric migration; T 1 , scaled time between split and migration event; T 2 , scaled time between migration event and present.
dispersal after the first year (Barthel et al., 2008). Given nest-site fidelity (Ridgway et al., 1991) in most populations, bidirectional movement of Smallmouth Bass cannot be ruled out. Regardless, individuals exhibiting the same movement strategies, either sedentary or migratory, are more likely to interbreed than individuals exhibiting different strategies (Barthel et al., 2008). These same life history traits need to be investigated in Neosho Bass to make a valid ecological comparison. Connectivity among populations may also be limited by impoundments. A tagging study in the Elk River basin showed that tagged Neosho Bass in Sycamore Creek, Elk River, and Buffalo Creek did not cross the reservoirriver interface created by Grand Lake O′ the Cherokee (Miller & Brewer, 2022). Alternatively, movement is likely limited by dams  in the Neosho range, creating bottlenecks and strong drift. the Elk River and the White River system, most likely owing to the karst topography. Several studies (Branson, 1963(Branson, , 1967 have found evidence of stream capture events in the CIH that may have occurred post-Pleistocene, a timeframe that could qualify as secondary contact with respect to the timing of species divergence. Intermittent periods of high water (flooding) may have also brought temporary stretches of connectivity. Admixture in the upper Arkansas River tributaries was best modeled by differential rates of asymmetric migration over two distinct epochs. While this demographic scenario has a distinct joint site frequency spectrum (Portik et al., 2017), we presume that the upper Arkansas River tributaries and Elk River populations have undergone similar natural processes given their geographic proximity. It is also likely that both populations have been subjected to anthropogenic introductions of Smallmouth Bass, either inadvertently or deliberately to promote angling (Johnson et al., 2009;Rahel, 2004). Varied timing and differential volumes of introduced fish in these two river systems could explain demography in these rivers. However, our inferences should be considered with caution, as we do not know of direct evidence of Smallmouth Bass stocking.
Admixture in both the Illinois River system and the Illinois Bayou and Big Piney Creek population was best explained by divergence with later secondary contact, more strongly implicating recent, anthropogenic introductions. We expected to find the Illinois River system has been subjected to secondary contact, because we know that it is admixed with the Skiatook Lake genomic cluster due to stocking of Lake Tenkiller with the same hatchery strain in the 1990s . We were surprised to obtain the same demographic history for the Illinois Bayou and Big Piney Creek population, because we inferred a greater amount of genetic drift following admixture in this part of the range from our admixture mapping results. Although supported gene flow models were consistent with widespread introductions of Smallmouth Bass from local sources, we have very little direct evidence of Smallmouth Bass stocking. Because we cannot determine the timing of secondary contact using the parameters derived in our demographic analysis, more data are needed before implicating recent introductions as the cause of admixture. Analysis of a full Smallmouth Bass genome and associated linkage map is needed to ascertain precise estimates of admixture timing.

| Conservation implications
Global rates of species loss, largely owing to anthropogenic habitat degradation and climate change, continue to climb at a scale warranting designation as Earth's sixth mass extinction (Ceballos et al., 2017;Kuipers et al., 2019). The projected consequences of these trends are dire, as biosphere stability depends on the interconnected ecological roles-e.g., productivity, predation, competi- Smallmouth Bass are keystone apex predators (Scott & Crossman, 1973) and obligate hosts for the larval stage of several freshwater mussels (Haag et al., 1999;Hoffman, 1999 It is urgent to predict the long-term fitness outcomes of gene flow in Neosho Bass and other taxa subjected to various forms of secondary contact, especially as species introductions continue to increase (Pearson et al., 2021). Gene flow may cause outbreeding depression through epistatic incompatibilities between derived alleles (Bateson, 1909;Dobzhansky, 1934;Muller, 1942) or undermine coadapted gene complexes that have evolved in isolation (Altukhov & Salmenkova, 1987;Goldberg et al., 2005;Moyle et al., 1986). Mixing may also have the opposite effect, facilitating heterosis by alleviating the genetic load of deleterious genes (Alleaume-Benharira et al., 2006), establishing stable tension zones (Arnold & Martin, 2010), or reversing stochastic loss of heterozygosity when adaptive alleles flow into small, genetically homogenous populations (Fitzpatrick et al., 2016;Hedrick, 1995;Tallmon et al., 2004;Willi et al., 2007). Even under the latter scenario of adaptive introgression, the ultimate result could be the loss of a distinct lineage through genetic swamping.
Neosho Bass are likely experiencing human-mediated hybridization and introgression due to introductions. Further research to determine whether introduced alleles are adaptive and becoming more prevalent in Neosho populations, potentially leading to genetic swamping, or, alternatively, if they are maladaptive and reducing relative fitness, would help evaluate the species' evolutionary trajectory. Highly complex patterns of diversification and gene flow have likely gone undetected in other species, both terrestrial and aquatic, that have evolved in variable environments and have been subjected to human-mediated introductions. This study offers a potential road map for conducting future analyses on nonmodel and potentially threatened species and will aid in the preservation of biodiversity.

ACK N OWLED G M ENTS
This project would not have been possible without the funding, resources, and support of the Missouri Department of Conservation.
We thank A. Miller, M. Thomas, J. Chilton, and Z. Morris for assistance with sample collection. We thank the staff at Floragenex, Inc., especially J. Walsh and J. Boone, for generating the genomic data analyzed in this study. We thank C. Chang, M. Malinski, M. Lipson, J.

FU N D I N G I N FO R M ATI O N
This project was funded by the Missouri Department of Conservation.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.